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Abstract 

This  paper  represents  orbit  propagation  and  detennination  of  Low  Eearth  Orbit(LEO)  satellites, 
which  is  immediate  orbit  satellite  positioning.  Beiseds,  satellite  global  positioning  system  (GPS) 
configured  satellite  receiver  provides  position  and  velocity  measures  by  navigating  the  filter  to  get 
the  coordinates  of  the  orbit  propagation  (OP),  while  the  satellite  orbit  coordinates  is  an  important 
basis  for  the  task  operation. 

The  main  contradictions  in  real  time  orbit  which  is  determined  by  the  problem  are:  orbit 
positioning  accuracy  and  the  amount  of  calculating  two  indicators;  how  to  balance  these  two 
indicators.  This  article  is  dedicated  to  solving  the  problem  of  tradeoffs.  To  plan  to  use  a  non-linear 
filtering  method  for  immediate  orbit  tasks  requires  more  precise  satellite  orbit  state  parameters  in  a 
short  time.  Although  the  traditional  extended  Kalman  filter  (Extended  Kalman  Filter,  EKF)  method 
is  widely  used,  its  linear  approximation  of  the  drawbacks  in  dealing  with  nonlinear  problems  was 
especially  evident,  without  compromising  Kalman  filter  (Unsented  Kalman  Filter,  UKF).  As  a  new 
non-linear  estimation  method,  it  is  measured  at  the  estimated  measurements  on  more  and  more 
applications.  This  paper  will  be  the  first  study  to  UKF  microsatellites  in  LEO  orbit  in  real  time, 
trying  to  explore  the  real-time  precision  orbit  determination  techniques. 

Through  the  preliminary  simulation  results,  they  show  that  based  on  orbit  mission  requirements 
and  conditions  using  UKF,  they  can  satisfy  the  positioning  accuracy  and  compute  two  indicators. 

Keyword:  global  positioning  system,  low  Earth  orbit  satellites,  orbit  propagation,  orbit 
determination,  Unsented  Kalman  Filter 

1.  Introduction 

The  satellite  orbit  detennination  (OD)  estimates  discrete  observating  the  position  and  velocity  of 
an  orbiting  object.  The  set  of  observations  includes  the  measurements  from  the  space  based  GPS 
receiver  (GPSR)  that  is  located  on  the  object  itself.  Satellite  orbit  propagation  (OP)  estimates  the 
future  state  of  motion  of  an  object,  whose  orbit  has  been  determined  from  past  observations.  The 
satellite’s  motion  is  described  by  a  set  of  approximate  equations  of  motion.  The  degree  of 
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approximation  depends  on  the  intended  use  of  orbital  infonnation.  Observations  are  subjects  to 
systematic  and  random  uncertainties;  therefore,  orbit  determination  and  propagation  are 
probabilistic. 

The  satellite  is  influenced  by  a  variety  of  external  forces,  including  terrestrial  gravity, 
atmospheric  drag,  multi-body  gravitation,  solar  radiation  pressure,  tides,  and  spacecraft  thrusters. 
Selection  of  forces  for  modeling  depends  on  the  accuracy  and  precision  required  from  the  OD 
process  and  the  amount  of  available  data.  The  complex  modeling  of  these  forces  results  in  a 
highly  non-linear  set  of  dynamical  equations.  Many  physical  and  computational  uncertainties 
limit  the  accuracy  and  precision  of  the  object  state  that  may  be  determined.  Similarly,  the 
observational  data  are  inherently  nonlinear  with  respect  to  the  state  of  motion  of  the  object  and 
some  influences  might  not  have  been  included  in  models  of  the  observation  of  the  state  of  motion. 

The  remainder  of  this  paper  is  organized  as  follows:  Section  2  describes  the  methodology  of 
GPSR  based  orbit  detennination;  Section  3  is  brief  introduction  of  the  disturbance  mathematical 
model;  Section  4  legends  the  orbit  determination  algorithm  description;  Section  5  describes  the 
GPS  measurement  models;  OP  algorithm  settings  description  in  the  Section  6;  and  Section  7  offers 
the  conclusion. 

2.  Methodology  of  GPSR  Based  Orbit  Determination 

Three  basic  strategies  are  presently  in  use  to  determine  precise  LEO  orbits  with  GPS.  They  are 
the  dynamic,  the  kinematic  or  non-dynamic,  and  the  hybrid  or  reduced-dynamic  strategies. 

The  dynamic  orbit  determination  approach  [1]  requires  precise  models  of  the  forces  acting  on 
user  object.  This  technique  has  been  applied  to  many  successful  space  vehicle  missions  and  has 
become  the  mainstream  of  Precision  OD  (POD)  approach.  Dynamic  model  errors  are  the  limiting 
factor  for  this  technique,  such  as  the  geopotential  model  errors  and  atmospheric  drag  model  errors, 
depending  on  the  dynamic  environment  of  the  user  space  vehicle.  With  the  continuous,  global, 
and  high  precision  GPS  tracking  data,  dynamic  model  parameters,  such  as  geopotential  parameters, 
can  be  tuned  effectively  to  reduce  the  effects  of  dynamic  model  error  in  the  context  of  dynamic 
approach.  The  dense  tracking  data  also  allows  for  the  frequent  estimation  of  empirical  parameters 
to  absorb  the  effects  of  unmodeled  or  mismodeled  dynamic  errors. 

The  kinematic  or  geometric  approach  does  not  require  the  description  of  the  dynamics  except 
for  possible  interpolation  between  solution  points  for  the  user  object,  and  the  orbit  solution  is 
referenced  to  the  phase  center  of  the  on-board  GPS  antenna  instead  of  the  space  vehicle’s  center  of 
mass.  Yunck  and  Wu  [2]  proposed  a  geometric  method  that  uses  the  continuous  record  of  object 
position  changes  obtained  from  the  GPS  carrier  phase  to  smooth  the  position  measurements  made 
with  pseudorange.  This  approach  assumes  the  accessibility  of  P-codes  at  both  the  LI  and  L2 
frequencies.  Byun  [3]  developed  a  kinematic  orbit  determination  algorithm  using  double-  and 
triple-differenced  GPS  carrier  phase  measurements.  Kinematic  solutions  are  more  sensitive  to 
geometrical  factors,  such  as  the  direction  of  the  GPS  satellites  and  the  GPS  orbit  accuracy,  and  they 
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require  the  resolution  of  phase  ambiguities. 

The  previous  two  strategies  each  have  counterbalancing  disadvantages:  various  mismodelling 
errors  in  dynamic  OD,  and  GPS  measurement  noise  in  kinematic  OD.  A  hybrid  dynamic  and 
kinematic  OD  strategy  would  down-weight  the  errors  caused  by  each  strategy,  but  still  utilize  the 
strengths  of  each.  One  such  strategy  has  been  devised  and  is  referred  to  as  reduced  dynamic  orbit 
determination.  The  reduced-dynamic  approach  [1]  uses  both  geometric  and  dynamic  information 
and  weighs  their  relative  strength  by  solving  for  local  geometric  position  corrections  using  a 
process  noise  model  to  absorb  dynamic  model  errors. 

2.1  Orbit  Propagation  Algorithm  Description 

The  orbit  propagation  algorithm  can  be  divided  into  two  main  tasks:  orbit  determination  and  orbit 
prediction  (propagation).  The  general  diagram  of  Orbit  Propagation  algorithm  is  described  on  the 
following  figure: 
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Fig.  1:  Orbit  Propagation  algorithm  diagram 

The  Orbit  Determination  algorithm  is  based  on  Unsented  Kalman  Filter  (UKF)  and  estimates  the 
object  state  vector  xA+1  and  covariance  matrix  Pt+1  from  discrete  observations.  The  set  of 
observations  include  the  measurements  (raw  data}  from  the  space  based  GPS  receiver  that  is 
located  on  the  space  vehicle  itself.  The  Orbit  Detennination  algorithm  includes  the  Orbit  Prediction 
task  as  Time  Update  stage  of  UFK. 

Orbit  Prediction  algorithm  calculates  the  future  state  of  motion  of  a  vehicle  xA+l  whose  orbit  has 
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been  detennined  from  past  observations.  Moreover  the  covariance  matrix  P/+l  is  propagated.  A 
numerical  integration  of  the  Dynamic  Model  is  applies  for  orbit  prediction. 

The  OP  solution  tk+l,xk+1  is  output  data  of  Orbit  Propagation  algorithm.  The  OP  solution  and 

covariance  matrix  can  be  obtained  as  from  Prediction  task  as  from  Determination  task.  The 
following  external  data  are  required  for  OP  solution  calculation: 

•  init  time  and  state  vector  tinit,xinit  for  algorithm  initialization  /  reinitialization; 

•  the  time  moment  tk+l  to  new  OP  solution  xi+1  calculation; 

•  the  set  of  observations  N  ^ ,  {raw  data}  for  new  estimation  x/;+1  calculation. 

The  following  input  data  are  obtained  from  the  previous  OP  solutions  calculation: 

•  the  last  OP  solution  tk ,  xk ; 

•  the  time  test  of  last  calculation  of  estimation  xi+1; 

•  the  covariance  matrix  P;  . 

The  maximal  time  of  continuous  propagation  T{^  ,  maximal  integration  time  step  hmax ,  minimal 
count  of  available  Navigation  SV  N  ^in ,  default  covariance  matrix  Vdef  are  used  for  algorithm 
control. 

2.2  Dynamic  Model 

A  dynamic  model  of  the  object  motion  essentially  adds  a  priori  knowledge  from  the  equations  of 
the  orbital  motion  to  the  kinematic  position  knowledge  as  obtained  from  the  raw  GPS  measurements. 
In  our  case,  the  dynamic  model  incorporates  the  complex  Earth  gravity  field  (EGM  96)  truncated  to 
order  and  degree  18.  Furthennore,  the  Sun  and  Moon  gravitation  and  atmospheric  drag  are 


accounted. 
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The  differential  dynamic  equation  of  motion  is  given  by: 
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where: 

vx,vy,v_-  are  the  ECEF  velocity  components  of  object; 

x,y,z- are  the  ECEF  radius  vector  components  of  object; 
b  -  is  the  receiver  clock  bias; 
d  -  is  the  receiver  clock  drift; 

a Nbody  —  is  Sun  and  Moon  gravitation  forces; 


(1) 


aGEo~  is  the  acceleration  due  to  geopotential ; 

Fdrag  =  {Fx,  Fy,  Fz }  -  is  a  perturbing  force  due  to  aerodynamic  drag, 

-  is  the  angular  velocity  of  Earth  rotation,  hire  and  belowQ^  7.2921151467e-5  rad/sec. 

The  user  coordinates  are  in  the  rotating  Earth- fixed  frame  (ECEF).  Although  this  adds  some 
complexity,  especially  due  to  the  Coriolis  and  centrifugal  acceleration  in  the  dynamic  model,  no 
reference  system  transfonnations  are  required  in  the  main  program  since  input  (initial  position  and 
velocity),  as  well  as  the  OP  output  are  consistently  referring  to  the  Earth-fixed  frame.  In  this  way, 
reference  system  transfonnations  may  completely  be  encapsulated  in  the  dynamic  model.  Moreover, 
some  dynamic  algorithms,  which  compute  the  accelerations  due  to  the  Earth’s  gravity  field  and  the 

atmospheric  drag,  may  be  formulated  simpler  in  an  Earth-fixed  than  in  an  inertial  frame. 

The  integration  is  canied  out  by  using  the  simple  fourth  order  Runge-Kutta  algorithm  without  any 

mechanism  of  step  adjustment  or  error  control.  The  fourth  order  Runge-Kutta  is  considered  an 
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adequate  numerical  integrator  due  to  its  simplicity,  fair  accuracy,  and  low  computational  burden. 
Numerical  integration  is  perfonned  in  the  rotating  Earth-Fixed  frame  (ECEF). 

Solar  radiation  can  be  neglected  because  its  effect  on  total  object  acceleration  is  much  smaller 
than  acceleration  due  to  perturbing  geopotential,  the  third  body  forces  from  the  Sun  and  the  Moon, 
and  atmospheric  drag.  According  to  [12]  magnitude  ratio  of  atmospheric  drag  and  solar  radiation 
for  average  size  spherical  objects  with  cx  =2.4  moving  along  the  circular  LEO 
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We  can  see  that  for  altitudes  less  than  600  km  solar  radiation  pressure  is  significantly  smaller 
than  atmospheric  drag.  Furthermore  atmospheric  drag  decreases  with  altitude  and  it  becomes 
negligible  for  altitudes  higher  700  km. 


3.  Disturbance  Mathematical  Model 
3.1  Earth  Gravity 

The  gravitational  potential  function  for  the  solid-body  mass  distribution  of  the  Earth  is  generally 
expressed  in  terms  of  a  spherical  hannonic  expansion,  referred  to  as  the  geopotential  in  the 
Earth- fixed  reference  frame  (ECEF).  The  gravitational  potential  function  U  is  defined  as  [9]: 


U  = 


GM 


N  n 


R\ 
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n=2  m= 0  \rj 
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where: 

U  —  Gravitational  potential  function  (m2 / s2  ); 

GM-  Earth’s  gravitational  constant,  hire  and  below  GM=  3986004. 418e8  m3 / s2  ; 
r  -  Distance  from  the  Earth’s  center  of  mass  (in); 

R  -  Semi-major  axis  of  the  WGS  84  Ellipsoid  ,  hire  and  below  R  =  3986004.4 18e8  m; 
n,m  -  Degree  and  order,  respectively; 
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(p  -  Geocentric  latitude; 

A  -  Geocentric  longitude; 

C nm 5 S nm  ~  Normalized  gravitational  coefficients,  it  is  defined  in  EGM-96  model  [9]; 

Pm(sin  (f>)-  Normalized  associated  Legendre  function; 

Pnm  (sin  </>)-  Associated  Legendre  function; 

P„  (sin  (/>)-  Legendre  polynomial. 


Pnmi  Sin^)  = 


(n  -  in).  (in  +  \  )k 


(n  +  m )! 


1/2 


Pnm  (Sin  <f)  -  (cos  (f) 

P„  (sin  <P)  = 


d" 


Pnm(  Sin^) 


'  \P„  (sin  ^)] 


d  (sin  <p)' 

1  dn  (■  i  j,  i  V’ 

- 7 - MSUI  1 

2 "n\  dismay 


In  Eq.  (2)  and  (3)  k  =  1  for  m  =  0  and  k  =  2  for  m  ^0. 

The  series  is  theoretically  valid  for  r>R. 

The  acceleration  due  to  geopotential  can  be  defined  as; 
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Variables  r,  (j),  i  are  related  with  object  ECEF  radius  vector  components  x,  y,  z  by: 


,  2  .  2  .  2  . 
r  =  -\  x  +  y  +  z  ; 
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According  to  Eq.(6) 
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Projections  of  Earth  gravitation  force  in  ECEF  are 
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The  following  recurrence  equation  can  be  applied  to  cos  mi  and  sin  mi  calculation; 
cos ((m  +  l)i)  =  cos  mi  cosi  -  sin  mi  sin  i 

(9) 

sin((m  +  l)i)  =  cos  mi  sin  i  +  sin  mi  cosi 

3.2  N-Body  Perturbation 

The  gravitational  perturbations  of  the  Sun,  Moon  and  other  planets  can  be  modeled  with  sufficient 
accuracy  using  point  mass  approximations.  In  the  geocentric  inertial  coordinate  system,  the  N-body 
accelerations  can  be  expressed  as  [12]: 
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where 

G  -  the  universal  gravitational  constant 

M;.  -  mass  of  the  i-th  perturbing  body  (Sun  or  Moon) 

rt  -  position  vector  of  the  i-th  perturbing  body  in  ECEF 

r  —r  -  position  vector  of  the  i-th  perturbing  body  with  respect  to  the  object  mass  center  in  ECEF, 
i  -  planet  index,  i  =  S  for  the  Sun  and  i  =  M  for  the  Moon. 


The  values  of  the  Sun  and  Moon  position  vectors  rt  can  be  obtained  from  the  following  equations: 


rs  =  nCr‘s 


v  —  dCre 


(11) 


where 

Q  -  is  a  transfer  matrix  from  current  Equatorial  Earth  Centered  Inertial  Frame  to  ECEF  defined  as: 


Q  = 


cos  Clet  sinQ/  0 
-sinQ/  cosQ/  0 
0  0  1 


(12) 


Qe—  is  angular  velocity  of  Earth; 

t  -  is  time  in  seconds  from  the  beginning  of  current  sidereal  day  as  defined  below:: 

t  =  0.997269566329084  •  GSMT;  (13) 


GSMT  -  is  Greenwich  Sidereal  Mean  Time,  see  [21]  fo  detail; 

C  -  is  a  transfer  matrix  from  Ecliptic  Earth  Centered  Inertial  Frame  of  J2000.0  to  current  Equatorial 
Earth  Centered  Inertial  Frame  as  defined  in  the  following  equation: 
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(14) 


£  -  is  the  mean  obliquity  of  the  ecliptic  as  defined  in  [21]; 

res  -  is  radius  vector  of  Sun  mass  center  in  Ecliptic  Earth  Centered  Inertial  Frame  of  J2000.0 
defined  as: 

/  n  \ 

COSAs 

rs=Ps  sin4  (15) 

l  0  J 

ps  -  is  mean  distance  between  Earth  and  Sun  mass  centers 

ps  =^t/(1.00014-0.01673cos/'  -  0.00014  *cos2/’)  (16) 

AU-  astronomical  unit,  hire  and  below  A  U=  1.49597870ell  m; 

V—  is  the  mean  anomaly  of  the  Sun,  see  [21]  and  [22]  for  detail; 

Xs  -  is  the  ecliptic  longitude  of  the  Sun  defined  as: 

Xs  =  Ls +1.9171 -sin/' +  0.02- sin2/'  + 0.0003 -sin3/'  (17) 

Ls  -  the  mean  longitude  of  the  Sun  as  defined  below: 

Ls  =180  + ( 100.46457166  +  35999.37244981 -T)  (18) 

T  -  is  Julian  centuries  from  J2000.0; 

reM  -  is  radius  vector  of  Moon  mass  center  in  Ecliptic  Earth  Centered  Inertial  Frame  of  J2000.0 
which  is  defined  as: 
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cos  XM  cos  P M  ^ 
sin/lM  cos  PM 
v  si  npM  , 


Pm 


is  the  distance  between  Earth  and  Moon  mass  centers  defined  as: 


Q  =  1  +  0 .0545  cos  /  +  0 .0100  *  cos(/  -  2D)  +  0 .0082  cos  2D  + 
0.0030 cos 2*/ +  0.0009  cos(/  +  2D)  +  0.0006  cos(/'-2D)  + 
0 .0004  cos(/  +  /  ’  -  2D)  +  0 .0003  cos(/  -  /  ’) 


(20) 


XM—  is  the  Moon  geocentric  longitude.  It  can  be  defined  as: 

XM  =  L  +  6.289sin /  -1.274  sin(/  -2D)  +  0 .658  sin  2D +  0.214  sin  2/- 
0.186  sin /’-0. 114  sin  2F  -  0 .059  sin  (2/  -  2D)  -  0 .057  sin(/  +  /'-2D)  + 
0.053  sin (/  + 2D) - 0.046  sin(/'- 2D) +  0.041  sin(/-/’)- 0.035  sinD- 
0.031  sin  (/  +  /’) -0.015  sin(2D  -  2D)  +  0.01 1  sin(/-4D) 


L  =  Q  +  F  -  is  the  mean  longitude  of  the  Moon; 

J3M  -  is  the  Moon  geocentric  latitude  as  defined  in  the  following  equation: 

PM  =5.128  sin  F  +  0.281  sin  (/  +  F)  -  0 .278  sin(F-/)- 

0.173  sin(F-  2D)  +  0.055  sin(F  +  2D-/)-  (22) 

0.046  sin  (/  +  F  -  2D)  +  0 .033  sin  (F  +  2D)  +  0.017  sin(2/  +  F) 


/,  /’  ,F,D,Q  -  are  fundamental  arguments  of  Moon  motion  theory,  they  are  defined  in  [2 1  ] . 
All  equations  of  this  item  are  given  according  to  [21]  and  [22], 


3.3  Atmospheric  Drag 

A  near-Earth  space  vehicle  of  arbitrary  shape  moving  with  some  velocity  v  in  an  atmosphere 


will  experience  drag  force.  The  drag  force  acceleration  can  be  modeled  as  [12]: 


L  drag 
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where 

p  -  the  atmospheric  density; 
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vr  -  the  object  velocity  vector  relative  to  the  atmosphere; 
in  -  mass  of  the  object; 

Cd  -  the  drag  coefficient  for  the  object; 

A  -  the  cross-sectional  area  of  the  main  body  perpendicular  to  vr . 

C  A 

The  parameter  ax  =  — —  is  sometimes  referred  to  as  the  ballistic  coefficient.  It  is  varied  during 

m 

an  orbital  motion  due  to  an  object  attitude  and  an  object  mass  center  evolutions  and  others  factors. 
The  middle  (typically)  value  of  a  ballistic  coefficient  is  used  because  this  factors  are  unknown  for  OP 

YYl 2 

algorithm.  Hire  and  below  crt  =  0.0 1  — . 

kg 

Different  empirical  dynamical  atmospheric  models  can  be  used  for  computing  the  atmospheric 
density.  These  include  the  Jacchia  71  [15],  Jacchia  77  [16],  the  Drag  Temperature  Model  (DTM)  [17], 
DTM-2000  [18],  MSIS-90  [19]  and  NRLMSISE-00  [20]  and  others.  The  density  computed  by 
using  any  of  these  models  could  be  in  error  anywhere  from  10%  to  over  200%  depending  on  solar 
activity.  A  deal  settings  are  used  for  aforementioned  atmospheric  models  computation.  For 
example  geomagnetic  activity  index,  daily  and  average  solar  flux  index  and  so  on.  They  are 
fluctuated  during  orbital  flight  and  must  be  monitored.  Sizeable  density  errors  can  be  acquired 
otherwise.  Furthermore  all  abovementioned  models  require  appreciable  computation  resources. 

According  to  aforesaid  in  the  current  project  local  atmosphere  density  model  [12]  is  employed.  It 
is  rough  density  model  relative  to  dynamical  models  but  this  model  is  very  easy  for  computation  and 
requires  no  settings  monitoring. 

Equations  for  density  calculation  are  the  following: 

h-hx 

p  =  pxe  H  ; 

H  =  70000.  +0.075(7z  -  hj;  (24) 

h  =  r  -  R 

where 

/Zj  -  the  reference  height,  /Zj  =  500000  m; 
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fccr 

px  -  the  atmospheric  density  on  reference  height,  px  =  2.e- 1 3  — =(- ; 

nr 

h  -  is  the  object  height; 

H  -  is  the  height  scale  of  the  unifonn  atmosphere. 

r  -  Distance  from  the  Earth’s  center  of  mass; 

R  -  Semi -major  axis  of  the  WGS  84  Ellipsoid, 

4.  Prediction  Algorithm  Description 

Orbit  Prediction  algorithm  calculates  the  future  state  of  motion  of  a  vehicle  whose  orbit 
has  been  determined  from  past  observations.  Moreover  the  covariance  matrix  is  propagated. 

To  construct  the  future  object  trajectory  the  Orbit  Prediction  algorithms  uses  the  dynamic  equation  of 
motion  given  in  section  2.1.  This  fundamental  equation  of  mechanics  provides  the  dynamic 
constraint  governing  the  orbit  solution.  The  true  acceleration  at  any  instant  depends  on  the  space 
vehicle  position  and  velocity  at  that  instant,  and  on  many  other  parameters  that  characterize  the 
forces  at  work.  The  predicted  trajectory  is  then  generated  by  integration  of  the  Eq.  (25): 

tk+h 

v(tk+h)=  J  v(t)dt  +  v(tk) 

h 

tk+h 

'■</;  +//I-  f  l(lkfc  + /1 1,) 

r,  <25> 

~  ,k+fl 
b(tk+h)=  J  d{t)dt  +  b(tk) 

d(tk+h)  =  d(tk) 

where 

h  -  is  the  integration  step,  it  is  limited  byhmax  (see  Fig.  1  for  detail); 
v(tk  +h)  =  {vx ,  v  ,  vz }  -  is  the  predicted  ECEF  velocity  vector  of  the  obj ect; 
r(tk  +h)  =  {x,  y,  z}—  is  the  predicted  ECEF  radius  vector  of  the  object; 
b(tk  +  h )  -  is  the  predicted  receiver  clock  bias; 

d  (tk  +  h)—  is  the  predicted  receiver  clock  drift; 
v(tk),r(tk),b(tk),d(tk)~  define  last  object  state. 
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The  object  state  derivatives  vector  is  defined  using  dynamic  motion  model  which  is 
described  in  section  2.1.  As  the  numerical  integrator  we  will  use  a  simple  fourth  order 
Runge-Kutta  algorithm  without  any  mechanism  of  step  adjustment  or  error  control. 
Numerical  integration  is  performed  in  an  ECEF  reference  frame. 

The  covariance  matrix  propagation  is  defined  below. 

The  differential  dynamic  equations  of  motion  are  given  by: 

x  =  f(x,t )  (26) 

where 

x  =  (vx ,  v  ,  v_ , x,  y, z, b , d)T  is  a  state  vector  that  includes  the  spacecraft  position  and  velocity  vectors, 
and  the  receiver  clock  bias  and  drift. 

The  propagation  of  xk+1  from  the  previous  state  for  covariance  matrix  propagation  is 
generated  by  the  following  reduced  equation: 

tk+h 

xk+1  =  fm  (xk  )  =  xk  +  J  /,.  (x,  t)dt  =  xk  +  hfr  (x,  t )  (27) 

where 

h  =  tk+l  -  tk  -  is  the  integration  step,  it  is  limited  by  hmax  (see  Fig.  1  for  detail), 
xk  -  is  the  state  vector  from  the  previous  step, 

f(x,t)-  is  the  reduced  dynamic  model  of  notion  witch  is  defined  below. 

A 

4 

E 
x 

y 

z 

b 
d 
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4.1  Orbit  Determination  Algorithm  Description 

The  Orbit  Determination  algorithm  applies  an  UKF  to  estimate  the  state  vector  which  comprises  8 
components: 

•  obj  ect  velocity  v  =  {vx,  v  ,  vz } 


•  object  position  r  =  {x,y,z} 

•  receiver  clock  bias  b 

•  receiver  clock  drift  d 


b  d 


The  diagram  of  Orbit  Determination  algorithm  describes  on  the  following. 


The  following  process  and  measurement  models  can  be  established: 


=/(**)+w* 

(29a) 

Zk=h{xk)  +  vk 

(29b) 

The  variables  in  the  above  equation  will  be  described 
xk  is  a  system  condition  vector  in  the  k  moment 

/(•)  is  unscented  system  model 

wk  is  a  dynamic  mixed  signal  in  the  k  moment 

Zk  is  a  measuring  dynamic  vector  in  the  k  moment 

/?(•)  is  a  unscented  system  measuring  model 

vk  is  dynamic  measuring  mixed  signal  in  the  k  moment 

The  measurement  vector  is  denoted  by  zk ,  the  process  noise  wk  and  the  measurement  noise 
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vk  are  assumed  to  be  zero-mean  white  noise.  The  process  noise  covariance  matrix  and  the 
measurement  noise  covariance  are  given  by  Qt  and  R  ,  respectively. 

The  system  error  covariance  matrix  (?.  is  as  follows. 


r.}  r\  [Qj  j  j 

E\wjwj\=< 

1  jf  |0  ,i*j 


The  measuring  error  covariance  matrix  R  is  as  follows. 


(30) 


E{v,vTJ}  =  \Rl  J  ’ 

1  Jl  lo  ,i*j 


Here,  wk  and  vk  are  independent  and  unrelated 


E{wiv]}  =  0,  V  i, 


(31) 


(32) 


4.2  Unscented  Kalman  Filter  Processing 
1)  Producing  Standard  Points  and  Calculation  Value 

Jo  =  x; 


Xi  =  *k  +(sj{n  +  £)Pk  ). 

f  y  J 

i  =  l,2,---,n 

(33) 

Zt+n=x;-(yl(n  +  e)Pk)j 

wm  -  8 

(«  +  £■) 

i  =  1, 2,  •  •  • ,  n 

K  =  K+{  i-«,2+^) 

(34) 

w"‘  -  wc  -  l 

'  ''  2  (n  +  s) 

The  parameter  is  a  scaling  parameter  defined  as 

6*2  =  s\  {n  +  e3 )  -  n  (35) 

The  positive  constants  C  ,  i=  1,2,3  are  used  as  parameter  of  the  method,  a  priori  and  a  posteriori 


estimates  of  the  state  are  denoted  by  xk  and  xk . 
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2)  Time  Updating: 

Condition  Predicted  Value 


(g~k)=f((xk)) 

Condition  Predicted  Average 

1=0 

Covariance  Matrix 

K  =  Z  w;  \ (f; )  -  *,-][(?; )  -  x; f  +  a 

i=0  L  J  L  J 

3)  Observation  Updating: 

Observation  Measurement  Predicted  Value 

Observation  Measurement  Predicted  Average 

i=0 

PX}  and  Pyy  update 

i=0 

i=0 

4)  Calculating  Kalman  Gain  Value 
K  =  P  P' 

5)  Updating  Estimated  Value  to  Measurement  Value 

6)  Updating  Condition  Error  Covariance  Matrix 

p  =P-KPKt 

1  k+ 1  U  ^k1  yy  k 

4.3  Unscented  Kalman  Filter  Flow  Chart 

The  flow  chart  of  the  UKF  is  as  shown  in  Fig. 2 


(36) 

(37) 

(38) 

(39) 

(40) 

(41) 

(42) 

(43) 

(44) 

(45) 
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Figure  2  The  algorithm  flow  chart  of  UKF 

The  time  update  phase  of  the  UKF  includes  the  propagation  of  state  vector  \~k+]  from  the  previous 

object  state  and  the  state  covariance  matrix Pk+X .  It  is  defined  in  item  section  2.1. 

The  subsequent  measurement  update  adjusts  the  state  vector  x“+l  components  and  state  covariance 
matrix  Pk+l  to  best  fit  the  GPS  pseudorange  and  pseudorange  rate  measurement  data. 

4.4  The  meassurement  update  phase 

The  measurement  residual  and  sensitivity  matrix  are  found  by  forming  the  computed  observation 
equation. 

The  model  for  a  GPS  pseudorange  measurement  is  given  by: 


19 


(46) 


where 


y  p  i,Xk  1 1  ’  4+1  )  Pk+ 1  +  4c+l  +  Pk 


k+ 1 


44+1  ^ji^GPS  U 1 1  )  (foPS yk+ 1)  (ZGPS  ZA+1  ) 


(47) 


is  the  geometric  range; 

5ci+i ,  yk+ j ,  zj.+1  are  the  positional  states  of  the  user  object  at  reception  time; 
xGPS  ,  yGPS ,  and  zGPS  are  the  positional  states  of  the  z-th  GPS  satellite  according  to  item, 

bk+ j  is  the  receiver  clock  offset; 


J3'k+1  accumulates  all  unmodeled  errors. 


Using  the  abovementioned  non-linear  measurement  equation  the  sensitivity  matrix  is  given  by  the 
Jacobian  matrix  of  partial  derivatives  of  non-linear  measurement  vector  with  respect  to  the  state 


vector  x: 


0  0  0 


( XGPS  Xk )  _  (y  GPS  y  k)  _  \ZGPS  Zk ) 


i  o 


pk 


pk 


pk 


(48) 


The  model  for  a  GPS  Doppler  measurement  is  given  by: 
Tvfe+n4+i)=(14+i  -U+i)- 


'k+\ 


-  r, 


k+ 1 


+  d  —  s‘ 

^Llk+ 1  ck+ 1 


A'+l 


r, 


k+ 1 


where 

va+i  ’  rL i  are  the  z-th  GPS  satellite  velocity  vector  and  radius  vector  according  to  item; 
vk+1  is  the  user  object  velocity; 
rM  is  the  coordinates  of  the  user  object  at  reception  time; 


(49) 


dk+l  is  the  receiver  clock  drift; 


s'k+1  accumulates  all  unmodeled  errors  of  the  Doppler  observation. 

Using  the  abovementioned  non-linear  measurement  equation  the  sensitivity  matrix  is  given  by  the 
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Jacobian  matrix  of  partial  derivatives  of  non-linear  measurement  vector  with  respect  to  the  state 
vector  x: 


»;(»*.,)= 


dD.  DD,  dD,  DD,  DD,  dD, 


where 


dD  _ 

{xgps  x) 

do  _ 

(fere- J;) 

SVy 

P 

dD  _ 

\ZGPS  Z) 

dv: 

P 

dD  _ 

-(vfs-v. 

dx 

dD  _ 

dy 

dD  _ 

fr-vj' 

dz 

i  (-^Gra  x) 


2\ 


,P 


P 


P 


P" 


dvx  dv  dvz  dx  dy  dz 


0  1 


^vGra_v  ^tyGPs  y ) 

(■ 


/> 


,vgps  |  i  [yGPs  y) 

'\v  vy\  „  3 

l  P  P 


-fr-vj 


{yGPs  yfccps  z ) 


pj 


-(vr-vj 


iXGPS  X\ZGPS  Z ) 


P 


| VGPS —y  ^GPS  y\ZGPS  Z ) 

f 


P 


_i VGPS —v  |  ^  iZGPS  Z ) 


,P 


P 


(50) 


(51) 


where 

x gps ’ y gps’ z gps  repesent rk+l  —  {xGPS,yGPS,zGPS} , 


vfVfVf*  repesent Vk+l  =  {(v, )Gra , (v y )‘Gra , (vz )GPS } ; 

x,y,z  repesent  rk+1  =  !-r  £+1  ’  y  k+l  ’  Zk+1 1  ’ 

Vv ,  vy ,  vz  repesent v‘+i  =  {(vv  )‘.+1 ,  (vv  )\+1 ,  (vz  )\+1 } ; 

D  represents  Di . 

If  both  pseudorange  and  Doppler  measurement  are  used  the  sensitivity  matrix  will  be  composed  of 
H  and  Hv  matrices  in  the  following  way 


Hlt .  = 


where  H  and  H  are  matrixes  size  of  [Nf’’,  x  8]  defined  as 


(52) 
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(53) 


Hp  = 


'Hi  - 

ii 

£cf 

l 

i _ 

nsf 

_ 

The  measurement  residuals,  or  innovations  sequence  is: 

yp(xk+l,tk+l) 

yv(*k+ i’tk+1)_ 


Ayk+ 1  = 


yP(xk+i^k+\) 

yv(xk+i’tk+i) 


(54) 


where 

y p  (xi+1 ,  tk+1 )  and  yv  {xk+l ,  tk+1 )  are  the  measured  pseudo  ranges  and  pseudo  range  rates  witch  are 
computed  according  to  section;  y  (xjt+1 ,  tk+1 )  and  v,  (3ct+1 ,  tk+l )  are  the  predicted  pseudo  ranges 

and  pseudo  range  rates  witch  are  computed  according  to  section. 

The  measurement  update  phase  uses  the  Kalman  equations  to  incorporate  the  information  given  by 
the  measurements  themselves,  and  obtains  improved  estimates  of  the  state  and  of  the  covariance: 

Kt„  =  LX.,  K.LX.,  +  L.  f 

xk+i=Xk+1+Kk+1Ayk+1  (55) 

L,  =(I~  Kkt,Hk„  )PM  (/  -  KtaHM  f  +  Klt,RktlKlt,T 
where  Rk+l  is  the  discrete  measurement  noise  covariance  which  is  basically  a  measurement  weight 
matrix. 

The  QR-decomposition  algorithm  is  applied  to  inverse  matrix  calculation.  The  general  idea  of  this 
algorithm  is  described  in  item. 


5.  GPS  Measurement  Models 

The  basic  measurement  types  that  will  be  employed  in  the  current  project  are  GPS  pseudorange 

and  Doppler  in  LI  frequency.  The  equation  of  the  pseudorange  in  LI  frequency  is  given  by: 

Pk  =  pk  + 1  k  +  c\dtGPS  (tk )  —  dtjj  (tk  )]  —  sk  (56) 

where  Pk  is  the  pseudorange  in  LI,  p‘k  is  the  geometric  range  to  the  7-th  satellite  at  the 
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observation  instant  tk  is  given  by 


Pk  *\j \X GPS  X^k  ))  +  ,y GPS  ))  +  ZGPS  Z(C  )) 


(57) 


I'k  is  the  ionospheric  delay,  c  is  the  vacuum  speed  of  light,  dtGPS  (tk )  is  the  GPS  satellite  clock 
offset,  dtu(tk )  is  the  receiver  clock  offset,  tk  is  the  observation  instant  in  GPS  time,  and  sk  is  a 
remnant  error  supposed  random  gaussian. 

The  numerically  controlled  oscillator  (NCO),  which  controls  the  carrier  tracking  loop,  provides  an 
indication  of  the  observed  frequency  shift  of  the  received  signal.  This  observed  frequency  differs 
from  the  nominal  LI  frequency  because  of  Doppler  shifts  produced  by  the  GPS  satellite  and  user 
motion,  as  well  as  the  frequency  error  or  drift  of  the  satellite  and  user  clocks.  The  equation  of  the 


Doppler  shift  in  LI  frequency  is  given  by: 


D'k=-\ 


v,  —  V 

-P - ^-LOSl 


L, 


(58) 


where  v‘k  is  the  z-th  GPS  satellite  velocity  at  the  observation  instant  tk ,  vu  is  the  receiver  velocity, 


LOS'k  is  the  line-of-sight  to  the  z-th  GPS  satellite  at  t, ,  and  Ll=1575.42  MHz  is  the  transmitted 


frequency. 

The  Doppler  can  be  converted  to  a  pseudorange  rate  observation  given  by  the  following: 


k  1 


(59) 


where/is  the  receiver  clock  drift  in  m/s;  and  elk  is  the  error  in  observation  in  m/s. 

Another  possible  GPS  measurement  model  is  a  linear  combination  of  GPS  LI  C/A  code  and 

carrier  phase.  Since  both  data  types  are  affected  by  systematic  ionospheric  errors  with  the  same 

magnitude  but  opposite  signs,  their  arithmetic  mean  is  free  of  ionospheric  errors.  This  approach,  as 

proposed  by  Yunck  in  1996  [1],  removes  the  dominant  systematic  error  source  for  raw  GPS  data, 

which  may  amount  to  10-20  m  [5]  at  low  elevations.  As  a  matter  of  fact,  the  resulting  so-called 

GRAPHIC  data  (Group  and  Phase  Ionospheric  Calibration)  provide  a  low-noise  biased  range  with  an 
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accuracy  of  half  the  C/A  code  noise.  A  drawback  of  using  the  GRAPHIC  data  type  originates  from 
the  employed  carrier  phase  data  which  introduce  range  biases  for  each  of  the  twelve  receiver 
channels.  As  consequence,  twelve  range  biases  have  to  be  adjusted  as  part  of  the  estimation  process 
which  significantly  complicates  the  orbit  determination  algorithms.  Finally,  it  has  to  be  noted,  that 
GPS  broadcast  ephemeris  errors  with  a  mean  standard  deviation  of  about  3  m  (3D  position)  and  1  m 
(User  Equivalent  Range  Error,  UERE)  are  still  present  in  real-time  applications  [11],  if  no  counter 
measures,  such  as  the  upload  of  precise  ephemeris,  are  taken. 


6.  OP  Algorithm  Settings 
6.1  Integration  Settings 

The  maximal  time  of  continuous  propagation  T^ax  =  2400  seconds  (it  is  specified  in  0). 

The  maximal  integration  time  step  hmax=  30  seconds.  It  was  defined  according  to  the  table  bellow 
which  describes  maximal  Runge-Kutta  method  errors  respectively  to  integration  step.  The  period  of 
dynamic  model  integration  is  one  turn. 


h,sec 

altitude,  km 

500 

1200 

dR,m 

dV,  m/sec 

dR,m 

dVm/sec 

1 

6.00E-07 

8.00E-08 

4.00E-07 

4.00E-08 

10 

0.009 

9.00E-06 

5.00E-03 

4.50E-06 

20 

0.16 

1.60E-04 

0.09 

8.50E-05 

30 

0.9 

9.00E-04 

0.5 

4.50E-04 

50 

9 

9.00E-03 

5 

4.50E-03 

100 

150 

0.15 

80 

0.075 

6.2  Dynamic  Model  Settings 

The  maximal  half-interval  of  multi  body  acceleration  fixing  r™ax  =  30  sec.  It  was  defined 
according  to  the  table  bellow  which  describes  integration  errors  respectively  to  the  half-interval.  The 
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period  of  dynamic  model  integration  is  one  turn. 


altitude, 

km 

half-interval  of  fixing 

10 

30 

60 

300 

dR,m 

dV,  m/sec 

dR,m 

dV,  m/sec 

dR,m 

dV,  m/sec 

dR,m 

dV,  m/sec 

500 

0.169 

0.000185 

0.551 

0.000602 

1.168 

0.00128 

5.806 

0.00635 

800 

0.189 

0.000193 

0.615 

0.000631 

1.304 

0.00134 

6.489 

0.00666 

1200 

0.216 

0.000204 

0.706 

0.000668 

1.497 

0.00142 

7.454 

0.00705 

The  fixed  multi-body  acceleration  components  are  available  on  time  interval  tfh  ±  r'"*x . 
Where  tfix  is  time  of  acceleration  fixing. 


6.3  Estimation  Settings 

The  minimal  count  of  available  Navigation  SV  N^in  =  2  (it  is  defined  by  test  results). 

The  discrete  state-noise  covariance  matrix  Q,  the  discrete  measurement  noise  covariance  R  and  the 
initial  covariance  matrix  of  P0  can  have  different  components  values  and  structure  for  special 
receiver  application.  According  to  the  requirements  0  in  current  protect  them  can  be  defined  for 
example  as  the  following: 


0.25  0  0 

0  0.25  0 

0  0  0.25 

0  0  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 


0  0  0 

0  0  0 

0  0  0 

25 e2  0  0 

0  25e2  0 

0  0  25 e2 

0  0  0 

0  0  0 


0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
le6  0 
0  le3 


(60) 
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2e-7 

0 

0 

0 

0 

0 

0 

0 

0 

2e  -  7 

0 

0 

0 

0 

0 

0 

0 

0 

2e  -  7 

0 

0 

0 

0 

0 

0  = 

0 

0 

0 

6e-2 

0 

0 

0 

0 

0 

0 

0 

0 

6e-2 

0 

0 

0 

0 

0 

0 

0 

0 

6e-2 

0 

0 

0 

0 

0 

0 

0 

0 

5e  - 1 

0 

0 

0 

0 

0 

0 

0 

0 

5e  -  3 

RPr  = 

55.9/77 

-  pseudo  range  measurement  dispersion 

777 

Rd  =  0.037 - 1- -  pseudo  range  rate  measurement  dispersion 

sec 


7.  Simulation  and  Analysis 

Using  the  Kalman  algorithm  to  estimate  orbit  propagation  and  detennination  in  this  section 
have  been  proposed  to  simulate,  to  validate  that  the  derivation  of  the  formula.  Simulation  results 
are  shown  in  Figure  2,  the  initial  conditions  are  selected  at  the  beginning  of  the  track  after 
leaving  a  balance  within  200  seconds  after  convergence.  The  simulation  results  as  expected. 

Simulated  conditions: 

Micro-satellite  altitude  500  km,  longitude  108°  and  latitude  35°  ,  sampling  time  4sec, 
On-modulator  magnitude  =  2,  satellite  attitude  motion  trajectory  is  shown  in  Figure  3.  The  UKF  of 
direct  observation  equation  is  used  in  the  simulation. 

In  Figure  3(a),  time  response  of  micro-satellite  measured,  estimated  and  difference  between 
measured  and  estimated.  Deliberately  made  a  real  track  star  with  an  initial  value  is  not  the  same 
kalman  filtering,  200  seconds  after  the  kalman  algorithm  convergeswithin  200sec.  In  Figure  3(b), 
time  response  of  micro-satellite  velocity  measured,  estimated  and  difference  between  measured  and 
estimated.  Same  as  Figure  3(a),  willfully  made  a  real  track  star  with  an  initial  value  is  not  the  same 
kalman  filtering,  200  seconds  after  the  kalman  algorithm  convergeswithin  200sec. 
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Time  response  of  actual  Eular  angle[((>0  0Q  \|/Q  ]=  [-22  36  45]deg 


Time  response  of  estimate  Eular  angle  [c|>0  0Q  v|/_  ]=  [-21.7064  42.2116  45.3727]deg 


Time  response  of  Error  of  Eular  angle 


Figure  9(a) 


Time  response  of  actual  body  rate  dt/d  [  4>0  0Q  \|/Q  ]=  [5  -3  -6]deg/sec 


Time  response  of  estimate  body  rate  dt/d  [  c[>0  0Q  \j/0  ]=  [-1.0651  -4.072  -7.6284]deg/sec 


Time  response  of  body  rate  error 


Figure  9(b) 


Simulation  results  show  that  when  the  micro-satellite  attitude  and  orbit  to  maintain  balance, 
satellite  orbital  position  and  velocity  in  the  estimated  value  and  the  measured  via  GPS  satellite 
computer  considerable  amount.  Initial  value  change  in  a  relatively  small  error,  the  maximum  error  of 
10  °.  If  the  number  of  GPS  by  the  change,  the  number  of  GPS  consists  of  four  changes  into  three  or 
less  Error  is  relatively  large  when  the  attitude  changes,  the  maximum  error  is  instantaneous  32  °.  In 
the  academic  theory  and  engineering  practice,  a  systematic  analysis  is  generally  considered 
consistent  with  the  conclusion  that  it  is  feasible. 
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